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Calculation of the even-odd energy difference in superfluid Fermi 
systems using the pseudopotential theory 
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PACS 74 . 20 . Fg - BCS theory and its development 
PACS 67 . 85 . Lm - Degenerate Fermi gases 

PACS 71 . 15 . Dx - Computational methodology (Brillouin zone sampling, iterative diagonalization, 
pseudopotential construction) 

Abstract - The pseudopotential theory is extended to the Bogoliubov-de Gennes equations to 
determine the excess energy when one atom is added to the trapped superfluid Fermi system with 
even number of atoms. Particular attention is paid to systems being at the Feshbach resonance 
point. The results for relatively small particle numbers are in harmony with the Monte Carlo 
calculations, but are also relevant for systems with larger particle numbers. Concerning the addi- 
tional one quasiparticle state we define and determine two new universal numbers to characterize 
its widths. 



The pseudopotential theory has proved to be an im- 
portant tool in different areas, as for instance solid state 
physics (sec, e.g., [I]), quantum chemistry (see, e.g., [2]), 
metal clusters (see, e.g., [3]). A general survey of the 
method regarding the interest of researchers in a wide 
range of applications is given in [3]. For a recent review 
about the approximations in electronic structure energy 
see ref. [5]. The idea was first formulated by Hellmann 
[6] and appeared also in Gombas' work [7]. The quantum 
mechanical foundation was started by Fenyes [8] (for the 
early history of the pseudopotential theory see [9|fT0] ) . 

In the quantum mechanical treatment of many par- 
ticle systems it is crucial to reduce the task from an 
-^Vtotai-particle problem to an N e ff particle one, where 
N e ff < Ntotai, when N to tai > 1. In atoms, molecules 
and metals, e.g., N e ff is the number of valence electrons. 
The pseudopotential method works by taking the advan- 
tage of the fact that the large negative potential energy 
felt by a valence electron is almost completely cancelled 
by its large positive kinetic energy, which is due to the os- 
cillations of its wave function inside the core. The method 
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deals with pseudo wave functions which are smooth and 
easier to handle. One expects that similar situation arises 
in any fermion system put in an external potential when 
the analog of the valence electron problem in a broad sense 
is investigated. Characteristic differences also show up in 
case of trapped Fermi gases as a consequence of the short- 
range nature of the interaction between the atoms and in 
particular when the superfluid state is realized, especially 
at the Feshbach resonance. 

The main goal of the present paper is to formulate and 
use the pseudopotential method for trapped Fermi gases 
at zero temperature. For the sake of concreteness the in- 
ternal degree of freedom of the atoms will be assumed to 
be twofold, one can have in mind for instance the 6 Li iso- 
topes as elements of a two component Fermi gas. The 
BCS side of the BCS-BEC crossover will be considered, 
where the interparticle interaction can be characterized 
by a negative s-wavc scattering length a. We will focus in 
particular on the Feshbach resonance point, where univer- 
sal properties set in, since \a\ — > oo there, and the range 
of the interaction can be neglected (for a review see [TTj). 

We start by recalling the Bogoliubov-de Gennes (BdG) 
equations, which can be considered as the generalization 
of the Hartree-Fock equations to include pair correlations 



p-1 



A. Csordas et al. 



and conceived as equations of a mean- field theory [T21[T3] ■ 
They are as follows: 

[H + Uint] u n + Av n = £„«„, 

A*u„ - [H + Uint] v n = e n v n , (1) 

or in short 

^4>n = £ n 0n (2) 

with <f> n = (it„, v n ). The operator H can be written as 



h 

H=- — V 2 + U ext (r)-^ 
2m 



(3) 



where m stands for the mass of the atom, \x for the 
chemical potential and U ex t for the confining potential, 
which will be assumed to be harmonic in the applications. 
Uinti?) denotes the diagonal part of the mean- field po- 
tential, while A(r) is the pair potential of the mean field 
theory. In a more general framework both J7» n t(r) and 
A(r) contain contributions from correlation effects, but in 
the local density approximation (LDA), adopted also here, 
they are local quantities. We are interested in a system of 
particle number N + 1 with N being even. In particular 
we want to calculate the change in the ground state energy 
when adding one particle to the A^-particle system. Since 
the theory adopted here is not working with fixed parti- 
cle numbers it is more precise to speak about even- and 
odd-number parity states [12], but for sake of simplicty 
the therminology above is often used. 

It is assumed in the following that the iV-body prob- 
lem is solved. It might mean the solution of the BdG 
equations in the BCS mean-field approach, in which case 
Ui n t(r) is generally put zero [TTlfM] and only A(r) should 
be solved self-consistently. It is argued that in the mean- 
field description this is a consistent choice [15]. At uni- 
tarity, however, there exists an extension of the den- 
sity functional theory, which determines also Ui nt (r) self- 
consistcntly [TBHTBj . In our general considerations there is 
no need to specify which case is adopted. We do not use 
the results explicitely, suppose only their existence in prin- 
ciple and do not restrict ourselves to such particle num- 
bers for which calculations have been realized. We avoid 
the details of the direct interaction C7j nt by fixing C7j„ t 
and A using the solution of the iV-body problem. In a 
more precise treatment one ought to solve the (N + 1)- 
particle problem self-consistently. (In a normal system, as 
the electrongas in an atom for instance, this amounts to 
neglecting the core polarization by the valence electron.) 
In our model pseudopotential at unitarity universal con- 
stants will carry all the informations concerning the N- 
atom problem. Therefore, we do not need to specify a 
rcgularisation procedure, as in refs. [T4lfT8]. 

It will be important in the following that eqs. ([1]) pro- 
vide positive and negative eigenvalues as well representing 
the charge conjugation symmetry [THU^] . More concretely 
if <j)( + } = (u, v) is an eigenfunction with eigenvalue e > 
then = (— v*,u*) is an eigenfunction with eigenvalue 



-e. The orthonormalization condition 



d 3 r [u* n u m + v*v m ] = S„ , 



(4) 



leads to the conclusion that eigenfunctions belonging to 
eigenvalues of different sign are orthogonal. 

Let us denote by E^ N ' the ground state energy of the 
N particle system, then, see, e.g., ref. [12] 



E (N+i) _ E (N) 



N is even. 



(5) 



For /j, the chemical potential of the A^-body system should 
be taken. The calculation of e m in as the smallest positive 
eigenvalue of eq. © raises severe technical difficulties for 
large N, since the corresponding eigenfunction becomes 
highly oscillating. To avoid it and to make the determina- 
tion of e m i„ manageable even in such a case we generalize 
the pseudopotential method to the BdG equations by in- 
troducing the modified equation 



n + u.„ 



(6) 



where <p n stands for the pscudo wave function </>„(r) 
(u n (r), v n (r)), and U ps is the nonlocal pseudopotential 



u t {v) 
v t (r) 



K*(r'),<(r')) 



+ (e n + et) 



;f^)(- t (r'),^(r')) 



(7) 



Here e t is positive by definition. </>t(r) = (u t (r), u t (r))'s 
are selected of those solutions of eq. ([1]) that are involved 
in building the functions C/, n t(r) and A(r). 

Two remarks are in order here. One has to include 
states of positive and negative eigenvalues to maintain 
the "charge conjugation" symmetry. The summation may 
extend only to a restricted set to avoid divergences, for 
instance to states 



{{u t \H + U mt \u t ) + (v t \H + U lnt \v t )] < 0. 



(8) 



This choice will be assumed in the following. One can eas- 
ily see that the solution of eq. ([6]) is not unique, any linear 
combination of and 4>\ ^ (satisfying the condition ([8])) 
can be added to <f) n without changing the eigenvalue e„. 
As a consequence the most general form of U ps can be 
written as 



(+) 



^X^i}. o) 



One can easily convince oneself that 



F 



(e„±e t )(^ 



(10) 



are quite general functions. This expression of the 
pseudopotential is the generalization of the Austin, Heine 
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Sham proposal |20j . A natural requirement might be that 
the pseudo wave function be as smooth as possible (for- 
mally, it has been proposed to minimize the kinetic energy 
|21l22j , see also [5] ) . An obvious condition is that the pseu- 
dopotcntial should cancel a large part of the external and 
the local mean-field potentials (see e.g., [Tll2"Il2"3"] ) . 

In this spirit we choose the generalized Austin type [TJ 
[20] form 

U H U A 

^ ps ^ ps 

ut* -uS* 



with 



U 



H 



U ps = 



£(u t (r) [U(r')-fx]ut(r') 



(11) 



+v* t (r) [[/(r')-MK(r') 



U ps = -^'(^(r)[C/(r')-M]<(r') 
t 

-<(r) [U(r') - »)v t (r') 



(12) Fig. 1: Energy difference E {N+1) - E (N) - fi in units of huj as 
a function of N. Both axes are logarithmic, 5 = 1.16. Straight 
line shows an N 1 ' 9 behavior. 



(13) 



The prime on ^ means that the summation extends to 
states specified above and 



U = U ( 



ext ~l~ Ui n t • 



(14) 



This expression of the pseudopotential corresponds to the 
choice 



F t (+) (r) = -([/(r)- M )U(r),t, t (r) 



i~\r) = -(^(r)-M)(-<(r),t4(r)Y 



(15) 



(16) 



We consider in the following a trapped gas at the Fes- 
hbach resonance, where universal properties character- 
ize the system. Let's define the Thomas-Fermi region 
(TFR), where the conditions for the Thomas-Fermi ap- 
proach arc fulfilled (here and in the following, of course, 
the generalized Thomas-Fermi theory including paircor- 
relations is meant supposing that the particle number 
is large enough 2-1, 25J). An essential contribution to 
Y,t ( u t( r X( r ') + v t (r)vt (r')) arises only in the TFR then 
and it can be taken zero outside this region. Furthermore, 
it can be replaced to a good approximation, by a Dirac- 
delta function 6(r — r') within the TFR if the number 
of particles is large enough, which in other words means 
that the nonlocality of the pseudopotential is disregarded. 
Concerning JT3]) it is decisive that in LDA ut and vt are 
plane waves and the condition (jSJ) selects wave numbers 
in pairs k, — k, that altogether leads to U A = 0. The 
pseudopotential becomes diagonal in this approximation 



M - U(t) 






U(t) - v 



5(r-r'), re TFR 

otherwise. 
(17) 



This pseudopotential can be regarded as a model one, 
since the approximation breaks down near the border of 
the TFR. The total hamiltonian of the model reads as 

+ U ps = I A 2«J , r, TFR, (18) 

and in the region r ^ TFR: 
/ h 2 v 2 , „ 

"~2^r + u ^t - fj- 



n + u ps = 



o 



2m 





■ U ext + (i 



(19) 

Here we used the fact that Ui nt = and A(r) =0 in LDA 
outside the TFR. 
At unitarity 



A(r)=6( f ,-U ext (r)), 



(20) 



where 6 is a universal constant (see, e.g., [11]). In our 
calculations 5 is an input. In contrast, S is obtained as an 
output from the methods used in refs. [141ITS] . 

We have carried out the calculation for a spherical sym- 
metric harmonic oscillator trap potential 



U ext (r) = ^muj 2 r 2 



fix 



(21) 



where x = r/R with R being the Thomas- Fermi ra- 
dius. First we have looked for a spherically symmetric 
solution (i.e., I = 0) to get the smallest eigenvalue of 
eqs. (|6]). (fT8|l . (fT9| . In principle one can solve the BdG 
equations for r 6 TFR and r ^ TFR by supposing that <j> n 
is finite and vanishes as r — > oo, and by carefully match- 
ing the wave functions at the border of TFR. But from the 
numerical point of view it is more controllable the method 
of expansion in a given basis set. In the numerics we ex- 
pand both components of <p n in the 3D harmonic oscillator 
basis 4>nim(^)- In the isotropic case I and m are good quan- 
tum numbers. Thus, the expansion is over the different n 
(radial quantum number) values. Matrix elements of the 
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kinetic energy and U ex t in this basis are known exacly. 
Numeric integration is applied for the matrix elements of 
A(r) and U ps (which are nonzero for r G TFR). After 
calculating the matrix elements of Q + U ps a simple di- 
agonalization yields the eigenvalue e„ and the expansion 
coefficients of u n and v n . Special attention is paid to the 
size of the truncated basis {</>„; m (r)|0 < n < n max }. (In 
the calculation we choose for n max such a value for which 
the classical turning point of 0n majc Zm(r) is at least 1.5 
times bigger than the Thomas-Fermi radius. This ensures 
that we have enough basis functions which "feel" both 
regions r G TFR and r TFR.) In figs. HHS] we had 
rimax = 100 and tested that choosing n ma x for a slightly 
bigger value (cf. =150) the calculated e n , u n and v n 
are practically the same as for n max = 100. 

In fig. [TJ the energy difference £ , ( JV+1 ) — ijW — /i as 
given in eq. © is plotted. It follows the N 1 / 9 law pre- 
dicted by Son [55] . Concerning the accuracy of the prefac- 
tor, it is remarkable that even for a relatively small number 
of particles as N = 20 our result lies between the Monte 
Carlo findings [271128] and the result of the superfluid den- 
sity functional calculation [16] . (Note that at this particle 
number one expects that the extra particle is in the s- 
state). In fig. [2] the components of the eigenfunction are 
drawn. They have nodes, the appearence of which can 
be traced back to the fact that the operator (fTS)) . (fTT)]) 
couples the bare (normal) states of positive and negative 
eigenvalues. 

The Thomas-Fermi approach breaks down in the vicin- 
ity of the TF-radius R, which circumstance may ques- 
tion the applicability of the pseudopotential (fT?]) there. 
One can show, however, that the width of the solution 
of eqs. ©, (HI]), dHJ) is S r oc R( d/R) 4 / 3 , where d is 
the oscillator length d = (mcu), while the width of 
the surface region (where gradient corrections are impor- 
tant) scales as oc R(d/R) 4 [25]. At large particle num- 
bers R 3> d is valid, which ensures that the error is small 
when extending the Thomas-Fermi density to the border. 
This makes possible to determine two new universal num- 
bers. In fig. |3] the functions w u = A u (hu>/ /i) +2 / 3 and 
w v = A v (huj/ [i) +2 / 3 are drawn together with the numer- 
ically determined widths of u and v, in units of R. The 
fitting provides A u = 1.2774 and A v = 1.14616. 

For the solution in case of I ^ one needs a refinement 
of the model (Then, of course, e m i n is to be understood 
for the given I). First of all we quote that at resonance 

U(v)-i^=Uu ext (v)- f i), (22) 

where £ is a universal constant, whose recent values are 
0.372(0.005) [29] and 0.376(5) [30]. The complete cancel- 
lation of U(r) — fx by the pseudopotential occurs now only 
in a more restricted region, namely in between the zeros 
< ri < r- 2 < R of V(r, £) defined as 

V{T,Q = Uu^T)-ri + ^£l. (23) 
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Fig. 2: The u (upper part) and the v (lower part) compo- 
nents of the pseudo wave function as a function of the di- 
mensionless variable of i = r/R. Parameters: [i — lOOTiaj, 
emm = 4.437 hu, S = 1.16. 



Outside this region, but within the Thomas-Fermi radius 
R the total potential is V (r, £) in the radial equation while 
outside R it is V(r, 1). Semiclassically in the region v\ < 
r < T2 the radial kinetic energy is positive. Note that 
in the semiclassical treatment the factor 1(1 + 1) in the 
centrifugal energy ought to be replaced by (I + 1/2) 2 . We 
keep the quantum-mechanical expresion, however, to be 
able to extrapolate the results to small Z-values. 

The solutions of the radial equations corresponding to 
the potential discussed above have interesting features. 
The minimal excitation energy scales as ~ /i -1 ~ iV -1 / 3 
as predicted by Son |26j , but the prefactor depends on the 
value of £ and is smaller than the result one gets from the 
estimated centrifugal energy at the Thomas-Fermi radius 
R. A detailed analysis shows that this effect is mainly 
due to the two component nature of the wave function. 
Furthermore, for larger I- values the excitation energy is 
no more proportional to 1(1 + 1). At about I pa l maX7 de- 
fined by the condition n = r-z, the energy curve merges 
into the excitation spectrum for 1 = 0, which might de- 
tect a breakdown of the mean-ficld-typc theory at such 
excitation energies. 

One can show that by replacing the wave function of 
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Fig. 3: The width of u (+) and v (x) belonging to the pseudo 
wave function for e min as a function of /i. The widths are 
measured in unit of R and fi in Tiui. The solid and dashed lines 
show the w u and w v curves (see text). 

the (N + l)-th atom by the pseudo wave function the 
wave function of the total system does not alter. This 
requirement was the starting point in |31j working within 
the Hartree-Fock theory (see also ref. [5] for the general- 
ization to the (iV + 2)-particle problem). The ground state 
of the iV-particle system can be written as 

iv& )=n^ +) io)-n^ )+ i°)' ( 24 ) 
t t 

where a^~^ + and + are the quasiparticle creation op- 
erators in negative and positive energy states, respectively, 
with a'^' the corresponding destruction operators. |0) 
stands for the vacuum of the atoms. 

The state of the (N + l)-body system can be given as a 
one-quasiparticle state 

a+|¥ >. (25) 

The corresponding state written in terms of the pseudo 
wave function reads as 

(«++E^ +) ^ +) +E c t ( " r ^ )+ )l*o), (26) 

which is equivalent to (f2"5)l as can be seen from the expres- 
sion of l^o)- Note that normalization constants have not 
been included in the discussion above. The coefficients in 
([26]) are 

= (e n ±e t )(4 ±] \^ n ). (27) 

In conclusion we have calculated the energy of an extra 
particle in trapped Fermi gases at unitarity by generaliz- 
ing the pscudopotcntial theory to the BdG-type equation. 
The background iV-particlc problem has been treated 
within the Thomas-Fermi theory. This approximation 
could be improved by including the Weizsacker-type cor- 
rection as determined in case of the unitary system |25| . 



Along with such an extension it would be appropriate to 
take into account also the nonlocality of the pseudopoten- 
tial, a feature lost when obtaining (fTTj) . All these are left 
for future work. 

Some final remarks are in order. In our treatment all 
about the neighbouring even-number state are compressed 
into two universal numbers 8 and £. The recent value of 
£ is somewhat smaller than previously used in different 
matching processes, which might influences the other pa- 
rameter values. This makes our choice to take the bare 
atomic mass in our calculation reasonable. Concerning S 
we have taken its mean-field value by the same reason and 
by the expectation that it alters only slightly when using 
different models. 

Applications of more elaborated forms of the pscudopo- 
tcntial to other systems along with details of the present 
work will be published elsewhere. 

The present work has been partially supported by the 
Hungarian Scientific Research Fund under Grant Nos. 
OTKA 77534/77629 and OTKA 75529. 
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